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ABSTRACT 

We derive the distribution of flux density of a compact source exhibiting strong diffractive scintil- 
lation. Our treatment accounts for arbitrary spectral averaging, spatially-extended source emission, 
and the possibility of intrinsic variability within the averaging time, as is typical for pulsars. We also 
derive the modulation index and present a technique for estimating the self-noise of the distribution, 
which can be used to identify amplitude variations on timescales shorter than the spectral accumu- 
lation time. Our results enable a for direct comparison with ultra-high resolution observations of 
pulsars, particularly single-pulse studies with Nyquist-limited resolution, and can be used to identify 
the spatial emission structure of individual pulses at a small fraction of the diffractive scale. 
Subject headings: Methods: data analysis - Methods: statistical - pulsars: general - Scattering - 
Techniques: high angular resolution 



1. INTRODUCTION 
1.1. Noise and Scintillation 

Scintillation of compact astrophysical radio sources has 
proven to be a remarkable probe of both the turbulent 
interstellar plasma and the extreme conditions of the 
emission regions of these sources. The present work is 
primarily motivated by observations of pulsars at me- 
ter and decimeter wavelengths, for which the AU-scale 
scattering resembles a stochastic 'lens' with a potential 
resolving power of less than a nanoarcsecond. However, 
our results apply to any compact object that exhibits 
diffractive scintillation. 

The emission and scattering physics involve the in- 
teraction and mixing of many random variables. At 
the heart is the amplitud e-modulated noise emitted by 
the source (|Ricketti 119751 ). Indeed, even the amplitude 
modulation is stochastic, with individual pulses showing 
complex broadband profile structures, from bursty mi- 
crostructure to s trong and often non-Gaussian pulse-to- 
pulse variations flRickett et al.l 119751 : 1 Cairns et al.1 120011 : 
iKramer et all 120021 ) . This emission then propagates 
through the dilute, turbulent plasma of the interstellar 
medium (ISM), where it undergoes dispersion and Fara- 
day rotation, in addition to being strongly scattered. Fi- 
nally, an observer samples the electric field in the pres- 
ence of additive background noise, which often dominates 
the signal. Furthermore, the statistics of this random 
process evolve as the Earth, pulsar, and ISM move rela- 
tive to each other. 

Since the size of the scattering disk (£) scales ap- 
proximately as the square of the observing wavelength 
(A), the angular resolution afforded by the scintillation 
is #diff ~ A/S oc A . Thus, longer wavelengths permit 
superior resolving power. However, at longer wavelength 
each scintillation element contains far fewer samples, be- 
cause of the relatively smaller timescale (Aid oc A -1 ) 
and bandwidth (A^ ocA~ 4 ) of the scintillation pattern. 
It therefore becomes subtle to decouple the scintillation 
from the source and background noise. 
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To address these difficulties, we present analytical ex- 
pressions for the probability density function (PDF) of 
intensity after arbitrary temporal averaging of the spec- 
tra. We describe the case when single-pulse spectra are 
formed and then are averaged in groups of N adjacent 
pulses; we exclusively use N to denote such temporal 
averaging. We account for the decorrelation of the scin- 
tillation pattern during this averaging and derive the re- 
sulting modification to the intensity PDF. We also ac- 
count for spatially-extended source emission and show 
that the leading order PDF correction is identical, up 
to scale, to that of the decorrelation of the scintilla- 
tion pattern within the averaging time. In addition, we 
demonstrate that spatially-extended emission affects the 
measured PDF of intensity, even in the limit of N = 1 
averaged spectra. This feature presents the remarkable 
opportunity of measuring the emitting region sizes for 
individual pulses, or of subclasses of pulses. 

The sources and variety of noise are diverse, stem- 
ming from the background noise and the amplitude- 
modulated noise emitted by the source. Measured quan- 
tities such as intensity then correspond to random vari- 
ables, and their stochastic variation introduces additional 
noise. Ap plied to the source, this noise is the familiar 
self-noise (IDickdll946t I Kulkarnilll98l IZmuidzinas1l200l 
IGwinn & Johnsonll201lD~ but the background and scintil- 
lation parameters have similar limitations. Such noise is 
generally unbiased and scales with intensity. Our deriva- 
tions fully account for all these types of noise. We also ac- 
count for artifacts arising from instrumental limitations, 
such as quantization of the analog signal and quadrature 
downconversion. These types of effects introduce addi- 
tional noise as well as bias and distortion. 

To supplement our calculated PDFs, we also analyze 
expressions based on the moments of the distributions. 
In particular, we extend the modulation index, a tradi- 
tional estimator of source emission geometry, to include 
the effects of self-noise and pulsar amplitude variations, 
and we analyze a technique that estimates self-noise. Our 
expressions refine these simple tools for intrinsic emission 
inference. 
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1.2. Strategy for Comparison with Observations 

An observation well-suited for comparison with our 
models has an observing bandwidth B 3> A^d, an obser- 
vation duration t \ JS 3> Aid, and a spectral accumulation 
time much longer than the pulse-broadening timescale. 
For each pulse period, the on-pulse and off-pulse regions 
yield estimates of the signal and noise intensities, respec- 
tively. These estimated parameters fully define a model 
of the intensity PDF (see , which is readily compared 
with a histogram of the measured data. 

The residual structure in this comparison contains in- 
formation about the spatial extent of the emission re- 
gion and the evolution of the scintillation pattern, as 
well as artifacts from instrumental limitations and errors 
in the estimated on-pulse and off-pulse intensities. The 
on-pulse and off-pulse averages once again parametrize a 
model for this residual structure (see Q, which is quan- 
tified by a single fitted parameter for amplitude. This 
fitted parameter yields the transverse size of the emis- 
sion region, if the scattering geometry is specified. 

1.3. Comparison with Previous Work 

Most past work has focused on the regime of many 
samples of the spectrum averaged together: N — > oo. 
For such dat a, the effects of nois e can be approximated 
as Gaussian. I Gwinn et al.l (|1998l) gave the expected dis- 
tribution of intensity for a strongly scintillating source, 
includ ing the effects of a s patially-extended emission re- 
gion. I Gwinn et al.l (j2000f ) extended these forms to ac- 
count for Gaussian self-noise at small intensity, and es- 
timated the contributions of averaging in time and fre- 
quency by deriving their effects on the modulatio n in- 
dex. In an unpublished manuscript, iCord es (2000) pre- 
sented expressions for the PDF of scintillation gain in 
terms of a Karhunen-Loeve expansion with the coeffi- 
cients defined by a Fredholm equation, and suggested 
approximate forms determined by x 2 distributio ns with 
the ap propriate number of degrees of freedom. iGwinnl 
pool gave the distribution of interferometric visibility 
for a pointlike or ex tended source in strong scintillation. 
I Gwinn et al.l ()2011l ) presented expressions for self-noise, 
suitable for large N. 

The present work focuses on the complementary regime 
of single spectral realizations, and averages of N such 
spectra. We connect with these previous works through 
the asympt o tic fo rms of our distributions. Similarly, 
Ivan Stratenl i|2()()!): analyzed statistics of Stokes parame- 
ters when N = 1 and showed, for instance, that the de- 
gree of polarization is undefined in this case. His results 
describe the high signal-to-noise regime with no effects of 
propagation or source emission structure, as appropriate 
for the high-frequency observations of gia nt pulses from 
the C rab pulsar that motivated his work (Hankins et al. 
l2003h . 

1.4. Outline of Paper 

In <|2j we briefly review the theoretical descriptions of 
pulsar emission and interstellar scintillation, and estab- 
lish our assumptions about each. Then, in we present 
an expression for the exact intensity PDF of a scintillat- 
ing point source, as well as useful approximations. We 
next apply these approximations in Sj4]to estimate modi- 
fications to the PDF arising from the decorrelation of the 



scintillation pattern within the averaging time, effects of 
finite source emission size, and observational constraints 
such as finite scintillation averages; we then use these 
expressions to derive the limiting resolution for emission 
size inference. In we derive the modulation index and 
analyze a technique that identifies self-noise. Finally, in 
we summarize our results and outline the prospects 
for pulsar observations. 

2. THEORETICAL BACKGROUND 

2.1. Pulsar Radio Emission 

Pulsar radio emission is thought to arise from coherent 
curvature radiation from a rclativistic electron-positron 
plasma, flowing outward along open field lines from the 
stellar surface to the light cylinder. Charge accelera- 
tion occurs in a 'gap' of depleted charge density where 
the force-free state cannot be maintained. Such sites 
initiate a pair- production cascade, which forms a sec- 
onda ry plasma (|Sturrocklll971l : iRuderman fc Sutherland! 
11975( 1. This secondary plasma generates the ra- 
dio emission, which is then beamed, ducted, and 
refracted as it propagates through the upper mag- 
netosphere dBarnard fc Aronsl [19861: lArons fc Barnard] 
119861: iLvutikov fc Parikhl 120001) . The gap may be lo- 
cated close to the polar ca p region (the 'polar gap') 
(jArons fc Scharlemannlll979T) . or along the boun dary of 
the open and closed field lines (th e 'slo t gap'1 J Aronsl 
119831 ) . An 'outer gap' ([Cheng et al.lll986l : lRomam1ll996| ). 
situated close to the light cylinder radius, is the favored 
site for much high-energy emission and may also con- 
tribute to the radio emission. Differing locations of the 
gap may be responsible for some of the emission variety 
between pulsars. 

The superposition of many independent radiators pro- 
duces white Gaussian noise. Thus, although the radio 
emission process is fundamentally coherent, by groups 
of many particles, the bulk emission may still be ef- 
fectively spatially and temporally incoherent. Intrinsic 
physical processes (and the pulsar's rotation) then mod- 
ulate the envelope of this random field, motivating an 
' amplitude- mo dulated noise' description of the emission 
(|Rickettlll975l ). Our results assume this type of emis- 
sion, confined to a region that is a small fraction of the 
diffractive scale r^ = A-D/(27r£) ~ 10 6 m, where D is the 
characteristic Earth-scatterer distance. 

2.2. Interstellar Scattering 

Density inhomogeneities in the dilute, turbulent 
plasma of the ISM scatter the radio emission, leading 
to multipath propagation. For most pulsars observed at 
decimeter wavelengths, the phases between pairs of paths 
differ by many turns and are therefore uncorrelated; this 
is the regime of 'strong' scattering. Equivalently, the 
diffractive scale is much smaller then the Fresnel scale, 
r F = VXD - 10° m. 

In general, the timescale for the evolution of the diffrac- 
tion pattern at the observer ranges from seconds to hours, 
and so is effectively static over the scattering timescale. 
In this case, the effect of scattering is to convolve the 
pulsar signal with a 'propagation kernel.' The charac- 
teristic width of the average propagation kernel is the 
pulse-broadening timescale to. Dispersion contributes an 
additional frequency-dependent delay that is readily re- 
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versed (|Hankinslll971[ ). The motivation and limitations 
of approximating scatt ering as a convoluti o n are given in 
much greater detail by iGwinn fc Johnson! <f201lh . 

Although this picture of strong scattering of a point 
source is quite general, a description of the effects of a 
spatially-extended emission region requires specification 
of the scattering geometry. We present explicit results 
for the case of thin-screen scattering, as is observation- 
ally motivated (jWilliamsonlll972tlKomesaroff et al.lll972t 
I Hill et ahl I2005D . We give additional reductions by as- 
suming a square-law phase structure function; the gener- 
alization to Kolmogorov or otherwise is straightforward 
and changes only the scale of the PDF modifications that 
we derive. 

2.3. Instrumental and Observational Requirements 

In practice, spectra must be formed over a finite ac- 
cumulation time £ a cc; the convolutio n representation of 
scatte ring requires t acc to (see IGwinn fc Johnsonl 
(|2011l ) for additional details). The unique nature of pul- 
sars and their generally short duty cycles allow a par- 
ticularly elegant limiting option: the formation of single- 
pulse spectra that contain all pulsed power. Our descrip- 
tions are exact in this limit, and appropriate whenever 
iacc ^ ^o- 

We assume that such spectra are formed and then 
averaged in groups of N adjacent pulses. We also as- 
sume that the averaging timescale is much shorter than 
the scintillation timescale - the 'snapshot image' of 
INaravan fc Goodman! (|1989D - which subsumes the much 
weaker condition that the scattering material is approx- 
imately static during the accumulation of a single spec- 
trum. 

Frequency averaging is analogous but inherits addi- 
tional information from non-stationary signals. The in- 
trinsic pulse shape incurs spectral mixing via its associ- 
ated frequency-domain convolution; the average is then 
of correlated intensities (see £|A.2.2[) . Thus, a descrip- 
tion of frequency averaging requires the specification of 
individual pulse profiles, whereas a description of tem- 
poral averaging requires only the phase-averaged source 
intensity for each pulse. 

Lastly, we assume that the data explore a large repre- 
sentation of the full ensemble of diffractive scintillation: 
(B/Avd) x (^obs/Aid) 3> 1. Comparison with observa- 
tions is most powerful if the stronger condition B 3> Afd 
is satisfied so that the intrinsic amplitudes of individual 
pulses can be estimated. 

2.4. Notation 

Despite the number and variety of random variables 
involved in the emission and scattering processes, much 
statistical homogeneity exists; our notation highlights 
this foundational uniformity. We use z x to denote a cir- 
cular complex Gaussian random variable with unit vari- 
ance, indexed by x. Some common examples of such vari- 
ates are a single realization of the scalar electric field or 
the complex gain from propagation. The PDF of z x is ex- 
ponential in its squared norm, P(zx) = ■K~ l e~\ z "\ 2 , with 
respect to the standard complex metric dRefz^Jdlmfz^]. 
This norm is implicit for all PDFs presented with respect 
to complex random variables. 

We also encounter many exponential random variables, 
usually in connection with a scintillation "gain," which 



are thus denoted G or G„ with P(G) = er G for G > 0. 
Note that we use P() to generically denote a PDF with 
respect to the given variables and parameters, as in the 
previous examples. 

3. DISTRIBUTION OF INTENSITY 

We now derive the expected PDF of intensity when the 
spectra are averaged over groups of N consecutive pulses. 
We assume that the scintillation pattern is frozen during 
each average, and that the source emission is pointlike. 
In 21 we relax these assumptions and discuss practical 
limitations. Our results apply to scalar electric fields that 
have been baseband shifted and coherently dedispersed. 

We first derive an exact expression for the PDF of in- 
tensity in i j3.ll However, for N > 1 this expression has 
a removable singularity at the origin and is analytically 
burdensome. These features motivate various alterna- 
tives, such as the i.i.d. approximation and the gamma 
approximation, which we derive in the following ^3.2. 1 1 
and £13.2.21 respectively. We present a traditional rep- 
resentation, the Gaussian approximation, in M3.2.3I for 
comparison. Both the i.i.d. and gamma approximations 
are simple and powerful, and we heavily utilize them for 
the remainder of the paper. 

3.1. The Exact Intensity PDF 

The pulsar emits white Gaussian noise ej, modulated 
by a time- varying envelope yis/i' Here, I s is a con- 
stant characteristic scale of the source intensity, e, is nor- 
malized to have unit variance, and the modulation /j is 
power-preserving. To account for pulse-to-pulse varia- 
tions in the phase-averaged flux density, we also include 
a dimensionless pulse amplitude factor Aj, where j in- 
dexes the spectrum (i.e. the pulse number). We require 
no assumptions about the timescale of variability of the 
envelope fi or the nature of the pulse-to-pulse variations 
Aj. Our treatment therefore accommodates arbitrary 
variability of the pulsar, such as the possibility of cor- 
related pulse-to-pulse variations, log-normal amplitude 
statistics, or nanosecond-scale bursts, which merely mod- 
ulates the Gaussian self- noise. 

The propagation then acts to convolve this intrinsic 
emission with a scattering kernel, as discussed in $2.'2{ 
The kernel similarly consists of a power-preserving enve- 
lope gi that modulates white Gaussian noise r\i of unit 
variance. The propagated signal is then superimposed 
with the background noise vTaA, where I n is the ampli- 
tude of the noise, and ft is white Gaussian noise of unit 
variance. We assume that the background noise level 
is effectively constant over the averaging time, although 
this restriction is easily relaxed. 

The observed scalar electric-field time series Xi and its 
Fourier-conjugate spectrum are thus given by 

x % = tJA~T s [(/e) * (gn)] t + (1) 

Xi = y/AjI B (/*?). (g * V)i + V^nft- 

Here, we denote Fourier conjugate variables with a 
tilde. In short, this equation describes the amplitude- 
modulated noise of the pulsar (y/T^fi^i) convolved with 
the stochastic propagation kernel (^ry-j) and added to 
background noise (^/T^fti). Conventionally, these terms 
are labeled by their noise-free limits: the pulse profile 
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\fi\ 2 , pulse- broadening function \gi\ 2 , gated signal I s , and 
background I n . 

Of course, e*, fji, and /3j are mutually independent (cir- 
cular complex Gaussian) white noise. We see that a sin- 
gle spectral sample Xi is of the form AjI s ZfZ g + ^/T^z-^, 
where Zf, z g , and z\, are each circular complex Gaus- 
sian random variables with unit variance. If the scintil- 
lation is held fixed (e.g. z g = const.), then the intensity 
\xi\ 2 is drawn from an exponential distribution with scale 
Ij = Ajl s \z g \ 2 + I n . The PDF of the average of N such 
intensities is then given by Eq. IA1I 



/ 



N 



P(I;N\z g )=Nj2 
i=i 



\ 



7N-2 



N 



(2) 



N 



(Mz s \ 2 y- 



9 t \N-2 



AT 

j=i n - ^) 

£=1 

Observe that when we designate AT as a distribution pa- 
rameter, we implicitly include the background J n and 
source amplitudes {A,/ s }. 

All that remains is to incorporate the distribution of 
the scintillation random variable z g . We assume that z g 
explores a representation of the full diffractive ensemble. 
We introduce a 'scintillation gain' parameter G = \z g \ 2 , 
so P(G) = e~ G . Integrating the scintillation ensemble 
through G then provides the PDF of intensity: 

/•oo -G 

P(I-N)=N dG - (3) 



(IsG) 



N 

><E 



(AjlsG + Q 



N-2 



NT 

' (AjJ 3 G + / n ) 



N 

n 



This equation gives the expected PDF of intensity for the 
average of TV pulses with amplitudes {I s Ai} and fixed 
scintillation, in the presence of white background noise 
with mean intensity I n , if the spectra explore a represen- 
tative ensemble of diffractive scintillation. Examples of 
this distribution for various N are shown in Figure [U 

Numerical evaluation of the PDF can be challenging 
because of the singularity at G = 0, although the in- 
tegrand is not divergent. To address this difficulty, we 
employ variable-precision arithmetic libraries; the addi- 
tional precision is essential for N > 10. 

3.2. Approximating the Intensity PDF 

Although we have derived the exact expression for the 
intensity PDF (Eq. [3]), we now derive several useful ap- 
proximations. These approximations resolve the singu- 
larity at G = 0, uncover the N — > oo behavior, and 
present a simplified framework for deriving and interpret- 
ing perturbations of the PDF arising from effects such as 
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Figure 1. PDF of intensity for different pulse averaging numbers 
JV = 1, 2, 100, and oo. The source and background intensities are 
unity, as are all the pulse amplitude factors Aj ; the dashed line de- 
notes the mean of 2.0. Because of the pulse amplitude degeneracy, 
the i.i.d. representation (Eq. |4jl must be used (and is exact). 

a finite size of the emission region, decorrelation of the 
scintillation pattern within the averaging time, and er- 
rors in model parameters. 

3.2.1. The i.i.d. Approximation 

We derive our first approximation by replacing each set 
of N pulse amplitude factors Aj by its average A = (Aj). 
This replacement simplifies the statistics because the in- 
tensities that are averaged under such assumptions are 
independent and identically distributed (i.i.d.); thus, we 
refer to this approximation scheme as the i.i.d. approxi- 
mation. 

We emphasize that the i.i.d. approximation does not 
ignore amplitude variability of the pulsar, but does re- 
duce the variability to a moving average. The difference 
is critical because, for example, the i.i.d. approximation 
is exact when N = 1. 

After substituting the local amplitude average A, the 
PDF of intensity during a fixed scintillation gain G (i.e. 
Eq. [2]) is given by an Erlang distribution (Eq. IA2I) : 



P(I;N\G) 



N 



N 



I 



N-l 



-g AJ a G + /„ 



(4) 



(N - 1)\ (AI S G + I n f 

As for Eq. [3l the scintillation ensemble can then 
be included via its exponential density: P(I; N) — 
f™dG e~ G P(I; N\G). 

3.2.2. The Gamma Approximation 

We now extend the i.i.d. approximation to partially ac- 
count for the effects of variation in the pulse amplitudes 
Aj within sets of N pulses. To achieve this extension, 
we replace N in Eq. U by an effective number of degrees 
of freedom iVeff, chosen such that the mean and variance 
match those of the exact PDF. Because N e g is not neces- 
sarily an integer, the factorial in the Erlang distribution 
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must be replaced by a gamma function; this modified 
form is the gamma distribution. The mean and variance 
of this gamma distribution are 



Mgamma 
2 



gamma 



AI S G + I a , 
(AI S G + J n ) 2 



(5) 



For the exact distribution, given by Eq. [51 the mean is 
simply the average of the marginal means, and the vari- 
ance is the average of the marginal variances, divided by 
N: 



/Wet = ^2 {AjI s G + I n ) , 
1 N 

^oxact = J^p ^2 ( A J^ G + 



(6) 



3=1 



Matching Eq. [5] and Eq. H gives the correspondence ap- 
propriate for the gamma distribution: 



(7) 



(Ef=ii + ^W4) : 



Ef=i(l + ^sG// n ) 



< TV. 



The scintillation ensemble is again included via integra- 
tion over G with its exponential weight. 

We refer to this extension of the i.i.d. approximation as 
the gamma approximation. The gamma approximation 
partially accounts for the pulse-to-pulse variability, and 
substantially improves the i.i.d. approximation for N 3> 
1. 

3.2.3. The Gaussian Approximation 

By the Central Limit Theorem, Eq. [4] approaches a 
Gaussian distribution as N — > oo, with mean / = 
AI S G + I n and variance SI 2 = (AI S G + I n f /N = I 2 /N. 
Accounting for the amplitude variations that are ignored 
by the i.i.d. approximation gives the generalized form, 



SI 2 = 



2(1 



j \2 SAL 



/N; i.e. pulse-to-pulse am- 



plitude variations contribute twice their variance to the 
source self-noise. These relationships both reflect the fa- 
miliar fact that self-noise is inve rsely propor tional to the 
number of averaged samples N (Dickg [l946f L 

A traditional approximation strategy replaces 
P(I;N\G) with this matched Gaussian distribution; we 
refer to this scheme as the Gaussian approximation. 



0.6 
0.4 
0.2 
0.0 
0.2 




T 



Exact 
Gaussian 
i.i.d. 
Gamma 



o 

a 
a 
< 



o 
u 
u 
W 



0.00 
-0.05 

0.01 
0.00 
-0.01 



1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 
N = 2 1 








1 1 1 II 1 1 II 1 1 


III 

N=10 ~ 






— 1 — j — 1 — 1 — 1 — 1 — 1 — — 1 — 1 — h 


— | — l — 1 — l — 
N = 50 - 


1 " 1 1 1 1 1 1 1 1 1 1 


I , , , " 



Intensity I 

Figure 2. Exact PDF of intensity for N = 2, 10, and 50, and 
the three approximations discussed in ^3.2 1 The source and back- 
ground intensities are unity, and the N pulse amplitude factors are 
uniformly distributed between and 2. 



Figure [2] demonstrates the relative quality of these 
three approximations in the presence of substantial am- 
plitude variations. The Gaussian approximation is infe- 
rior to both the i.i.d. and gamma approximations, partic- 
ularly for small N; indeed, the Gaussian approximation 
includes nonzero probability for negative intensities. The 
i.i.d. and gamma approximations, on the other hand, give 
excellent representations for arbitrary N and are much 
easier to implement numerically than the exact represen- 
tation of Eq. [H particularly for large N. 

Our principal use of these approximations is to derive 
analytical estimates for PDF modifications arising from 
effects such as finite emission size and decorrelation of 
the scintillation pattern within the averaging time. We 
favor the relatively simple i.i.d. approximation for this 
purpose. While the small errors in the approximations 
might be significant in some measured PDFs, we expect 
the similarly small errors on the expected <1% PDF 
modifications to be negligible. Furthermore, the most 
intriguing degree of averaging is N = 1, which allows 
single-pulse studies and complete decoupling of emission 
size from decorrelation; both the i.i.d. and gamma ap- 
proximations are, of course, exact for this case. 

We therefore recommend the exact representation for 
PDF comparisons (or the gamma approximation if its ac- 
curacy suffices), and now derive PDF modifications using 
the i.i.d. approximation. 



3.3. Utility of Approximations 

These three approximations are particularly useful be- 
cause they reduce the N + 1 degrees of freedom in a 
single averaged spectra to three: an overall scale parti- 
tioned into source and background contributions, and a 
measure of the averaging. In fact, the i.i.d. approxima- 
tion has only two degrees of freedom because the degree 
of averaging is TV; this reduction comes at the expense of 
ignoring pulse amplitude variance. 



4. PDF MODIFICATIONS 

We now analyze effects that modify the PDF of inten- 
sity relative to the form derived in f}3] Modifications can 
arise from physical processes both intrinsic and extrinsic 
to the source, as well as from instrumental effects. Since 
the model PDF of intensity of <J3] depends only on direct 
observables (the average on-pulse and off-pulse intensi- 
ties), the residual structure in a measured distribution 
relative to the corresponding model reflects these modi- 
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fications. We also analyze the consequences of errors in 
the model parameters. 

We find degeneracies among these many types of ef- 
fects. Many of these degeneracies can be resolved by 
varying the degree of averaging N or the bandwidth used 
to construct the PDFs. We quantify the significance of 
each effect in terms of its expected bias on the inferred 
emission size 7 S and give strategies for identifying each 
source of bias. 

However, despite the nearly identical form result- 
ing from these varied effects, the derivation for each 
is unique. First, we account for decorrelation of the 
scintillation pattern within the averaging time in §4.11 
Then, we de rive the effects of spatially-extended source 
emission in ^4.21 These two effects are degenerate, to 
excellent approximation, and we analyze their relative 
strength in §4.31 Next, we calculate the consequences of 
both noise and bias in the model parameters in £14.41 and 
we estimate the influence of instrumental distortion in 
£14.51 Finally, we derive the expected spatial resolution 
afforded by these scintillation statistics in §4.61 



4.1. Effects of Temporal Decorrelation 

The scintillation pattern slowly evolves in response to 
the relative motions of the Earth, pulsar, and scattering 
medium, thereby affecting statistics of averaged intensi- 
ties. We now derive the effects of this evolution on the 
observed distribution of intensity, in terms of the tempo- 
ral autocorrelation function of the observed scintillation 
pattern. 

Although the scintillation pattern evolves somewhat 
similarly with frequency, the effects of frequency aver- 
aging fundamentally differ from those of temporal av- 
eraging. Namely, for temporal averaging, the intensi- 
ties within each average are independent, with correlated 
scales (i.e. the respective scintillation gains); whereas for 
frequency averaging, both the intensities within each av- 
erage and their scales can be correlated because of the 
intrinsic amplitude modulation. As noted in §2.31 spec- 
tral statistics after frequency averaging depend on in- 
dividual pulse profiles, whereas spectral statistics after 
temporal averaging only depend on the phase-averaged 
flux densities. 

Procedurally, we must index the scintillation factors G 
by pulse number. The resulting set {Gi] is then a col- 
lection of correlated exponential random variables, with 
covariance determined by the autocorrelation function of 
the measured dynamic spectrum: 



(GiGj) - 1 

{[I(t,u)-I n ] [I{t + M l3 ,v)-I n ]) 
(I(t,v)-I n ) 2 



(8) 



Details of this type of multivariate exponential distribu- 
tion are given in Appendix IA.2I 

In the 'snapshot image' regime, the averaging timescale 
is much shorter than the decorrelation timescale of the 
scintillation pattern, so ss 1. The leading order decor- 
relation effects in (1 — Fy) then provide an excellent ap- 



proximation. Using Eq. IA51 we obtain 
P(J;JV s {r tf })-P(J;JV) 



(9) 



- J2 (1 - r «) / d G Ge 



i<3 



d 2 P(I;N\{G t }) 
dG l dG ] 



Note that the partial derivatives are easily evaluated us- 
ing the characteristic function of P(I; N\{Gi}): ip(k) — 

\\ t (l — ilik/N^j , where Ig = AjI s Gj + I n . For exam- 
ple, applying the i.i.d. approximation (A, = A) yields 
the following relationship: 



d 2 P(I;N\{Gi}) 
dG l dG j 



1 



d 2 P(I;N\G) 



1) 



dG 2 



(10) 



G N(N- 

We can also reduce the sum over in Eq. |H] by 
approximating the evolution of the scintillation pattern 
as relative motion of the pulsar and observer with re- 
spect to a 'frozen' random screen. The decorrelation 
in time is then simply expressed in terms of the phase 
struc t ure function of the s cattering medium (jGoodmanl 
[19851: iCornwell et alJ [l989V For example, a square-law 

phase structure function has = exp — (Aty / Aid) 2 , 

even in the presence of axial anisotropy (jRickettl Il977t 
Uohnson et al.l l2012b[ ). If the averaging timescale 
is much less than the decorrelation timescale, then 
E,-<j (1 " r y ) « N 2 (N 2 - 1) / (l2Ar d 2 ), where Ar d is 
the decorrelation timescale in pulse periods. The quartic 
dependence on N leads to a rapid onset of decorrelation 
artifacts in the observed PDF after only modest averag- 
ing. 

Combining this approximation for the decorrelation 
behavior with the general form for the PDF modifica- 
tion (Eq. [9]) , and applying the reduction afforded by the 
i.i.d. approximation (Eq. 1 10[) then provides the following 
estimate: 



P(I;N, {Tij}) - P(I; N) 
N(N-1 
~ 12At, 2 



(11) 



dG Ge 



_ G d 2 P(P,N\G) 
dG 2 



4.2. Effects of Extended Source Emission 

A spatially-extended emission region smoothes the 
fluctuations from scintillation by superimposing many 
slightly offset copies of the diffraction pattern at the 
observer. In fact, extended source emission is merely 
one example in a broad class of physical effects that al- 
ter the distribution of scintillation gain. Other possi- 
bilities arise from modified scattering assumptions, such 
as weak scintillation or L evy-type scattering statistics 
(Bold vrev fc Gw inn 2003). To account for these cases, 
we introduce a generalized scintillation gain random vari- 
able G, which is no longer constrained to follow expo- 
nential statistics. Specification of P(G), combined with 
P(I;N\Q) given by Eq. [3l then provides the intensity 
P DF P(I;N) 

iGwinn et all (|1998ft examined P(Q) for thin -screen 
scattering of a Gaussian distribution of source intensity. 
They demonstrated that, if the source emission region 
is small relative to the magnified diffractive scale, Q is 
well- approximated by the convolution of three indepen- 
dent exponential random variables. The scales of these 
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random varia bles reflect the size of the emission region. 
IGwinnl (|2001[) derived an analogous result for interfero- 
metric visibility. 

By extending the description of §3.1[ we now demon- 
strate that an equivalent form arises for any spatially- 
incoherent, compact emission region in strong scattering. 
We parametrize the transverse coordinates of the source 
by s. Each emitting location has independent source 
noise, and the emission envelope at each location may 
also vary. However, the propagation kernel is correlated 
over the emission region. A single spectral sample thus 
takes the form 



d 2 s v /A(s)I s z f (s)z g (s) 1 + y/hz h . (12 



Xi — 



We specify our source coordinates, s, to be centered on 
the spatial mean of source intensity: J d 2 s sA(s) = 0. 

If the emission is confined to a small fraction of the 
magnified diffractive scale, then the scintillation random 
variable, z g (s), will be strongly correlated over the re- 
gion of source intensity. We therefore expand to linear 
order: zJs) sa z g (0) + (s • V)z g J s=0 . The source term 
in Eq. [T2] is then a convolution of three complex Gaus- 
sian random variables. Because of our choice of origin 
for s, these three random variables are mutually uncor- 
related during a fixed scintillation pattern. The marginal 
variances therefore add to give Q (up to overall normal- 
ization). In addition, the scintillation random variable, 
z g (s), is uncorrelated with its spatial derivatives, so the 
scales of the three variances are mutually independent. 
Thus, Q is the convolution of three independent expo- 
nential random variables: 



Q = 



l 



1 + 7s,l + 7s,2 



(G + 7s ,iGi+7 s ,2G 2 ). (13) 



The dimensionless subsidiary scales 7 S ,,: <C 1 contain 
information about the transverse extent of the source 
emission. The scaling prefactor is merely chosen so that 
(0) = 1 and could equally well be absorbed into the def- 
inition of source inte nsity. 

Using the results of lGwinn et al.l (|1998h , we now explic- 
itly relate the parameters 7^ to the emission geometry 
by assuming a Gaussian distribution of source intensity 
and thin-screen scattering. In this case, the 7 s ,i give the 
squared size of the source in orthogonal directions ii, in 
units of the magnified diffractive scale: 



D 



7s. 



R -J-A 

2tt 9i , 



(14) 



Here, D is the characteristic observer-scatterer distance, 
R is the characteristic source-scatterer distance, A is the 
observing wavelength, 9i is the angular size of the scatter- 
ing disk along ii , and <ii is the standard deviation of the 
distribution of source intensity along Xi. This represen- 
tation accommodates anisotropy of both the source emis- 
sion and the scattering. If the source emission arises from 
a circular Gaussian intensity profile, then the FWHM of 
the emission region is 2%/ln 4er 2.35c, where a = (Tj. 

Efforts to infer emission size often attempt to isolate 
the distribution P{Q) through sufficient spectral aver- 
aging; the modulation index m then qu antifies the ef- 
fects of emission size on this distribution (Sa rpeterlll967l : 



ICohen et al1fl967l: IGwinn et al.lll~998l) : 

_2 (S 2N 



(G) 2 1 

1 + 7s,l + 7s 2 2 

(1 + 7s. 1 + 7s, 2) 



(15) 



l-2( 7M +7 Si2 ). 



In g57TJ we discuss the modulation index in depth and 
extend Eq.[l5]to account for the effects of self- noise, tem- 
poral decorrelation within the averaging time, and pulsar 
amplitude variability. 

To determine the PDF modification resulting from a 
finite emission size, we first expand the PDF for a fixed 
scintillation pattern to linear order in j s y. 

P(I;N, ls \g)-P(I;N, ls =0\G) (16) 

dP{I;N, % =0\G) 



[Gi%,i + G 2 7s, 2 - G (7 S , 1 + 7s i2 )] 



dG 



Of course, the scintillation ensemble must now be in- 
cluded. Integrating each Gj with its respective weight 
e -Gi gj ves 



P(I;N, ls \G)-P(I;N, ls =0\G) 

dP(I;N, ls 



(17) 



2 7 s (1 - G) 



=0|G) 



dG 



Here, we have eliminated the degeneracy that appears 
at linear order between 7 S) j and 7 Sj2 by setting them 
both equal to a single characteristic size 7 S . If the effects 
of source emission size are substantial, then additional 
emission information, such as elongation or core/cone ge- 
ometry, can be inferred from the higher order terms, al- 
though higher terms in the original gain expansion must 
also be retained. To facilitate a comparison with the 
decorrelation effects derived in H4.ll we integrate by parts 
after combining with the final scintillation weight: 



P{I\ N, 7 S ) — P(I; N, 7 S =0) 



(18) 



-2 7s / dG Ge 

IQ 



_ G d 2 P(I;N, % =0\G) 
8G 2 



See Figure [3] for a comparison of these effects of emission 
size for various N. 

The approximate form of Eq. [T5] must be applied with 
some care; it fails in the limit N — ¥ 00, for example, be- 
cause of the resulting divergent derivatives of the delta 
function that characterizes the measured intensity. How- 
ever, the large- N limit is still easily accessible numer- 
ically by applying the full distribution of scint illation 
gam, as given by Eq. 30 in [G winn "etall GUI). Sub- 
stituting 7s = 7s, 1 = 7s, 2 in their expression and scaling 
to normalize the mean gives 



P{Q) 



l + 2 7s 
(l~7s) 2 



-(l+27s)0 



(19) 



^ + (l + S)-2g 7 s)e-( 2+ ^) e 
7s / 



When combined with the conditional PDF P(I;N\Q) 
(Eq. [2]) , this equation produces the scintillation averaged 
PDF of intensity. The disadvantage of this method is 
that the PDF must be calculated separately for every 
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Figure 3. PDF of intensity and the corresponding modification 
expected for an extended emission region, for N = l, 10, and 100. 
The source and background intensities are unity, as are all pulse 
amplitudes Aj. 

value of 7s. Fits for emission size effects are then non- 
linear and computationally expensive, in contrast with 
the simple linear fit enabled by Eq. [T8] Fortunately, the 
linear expansion of Eq. [18] is easily adequate for the full 
expected range of N, even exceeding 10 6 , if 7 S <C 1. For 
7s > 0.1, the linear approximation begins to differ sub- 
stantially from the exact expression, but this region re- 
quires the appropriate quadratic term added to the orig- 
inal gain expansion as well. 

Comparison of Eq. [18] with Eq. QTJ shows that, at lin- 
ear order and to the accuracy of the i.i.d. approxima- 
tion for the pulse amplitudes, the effects of an extended 
emission region on the PDF of intensity are identical to 
the effects of temporal decorrelation, up to overall scale; 
this relationship is further analyzed in H4.3I Moreover, 
the modification arising from extended emission increases 
with the number of samples averaged, as was the case for 
temporal decorrelation (see Figure [3]). These effects also 
depend on the gated signal-to-noise S = (A) I s /I n of the 
observation; see Figured] Indeed, the dependence on S, 
when combined with the Poisson noise, determines the 
resolution limit of an observation. We derive this limit 
explicitly in £14.61 

4.2.1. Effects of a Non-Static Emission Region 

Our analysis is well-suited for long spectral accumula- 
tion times, over which the pulsed emission may be highly 
non-stationary, and the pulsar motion may be substan- 
tial. For instance, the emission might be bursty; the 
received radiation would then represent a superposition 
of spatially-offset emission sites. The inferred emission 
size, as quantified by 7 S , will then be a characteristic size 
for the emitting region over the corresponding retarded 
time. This characteristic size weights the spatial stan- 
dard deviation by the integrated flux density at each lo- 
cation. Thus, spatially-offset emitting regions contribute 



to the characteristic size, if they emit within the same 
accumulation time. In this sense, our estimate presents 
an upper limit on emission size, since the emitting region 
at a single instant may be small relative to the overall 
emission site. Furthermore, a small emitting region may 
undergo substantial displacement over the accumulation 
time, from lateral motion of the pulsar or rotation; both 
of these effects will present an upward bias of the size 
estimate. However, because of the weighting by flux den- 
sity, weak extended emission may be overwhelmed by a 
strong pointlike component, yielding a nearly zero char- 
acteristic size. Any interpretation of a measured 7 S must 
address these competing factors. 

4.3. Comparison of Emission Size and Temporal 
Decorrelation 

As we have already observed, comparison of Eq. [TT] 
with Eq. [18] demonstrates that the linear corrections 
to the intensity PDF for effects of finite emission size 
and temporal decorrelation are perturbations of identi- 
cal shape, and differing weights, at least within the ex- 
cellent i.i.d. approximation for the averaged intensities. 
The relative strength of these two perturbation weights 
is 



size 



decorrelation 



24 



N(N- 1) 



7s 



(20) 



where we have substituted the square-law autocorrela- 
tion given in £14.11 Eq. [2UJ allows immediate assessment 
of the relative importance of effects. For example, if the 
source extends over 5% of the diffractive scale, then ef- 
fects arising from the finite emission size will dominate 
those from decorrelation of the scintillation pattern as 
long as the averaging is over less than a quarter of the 
diffractive timescale. 

4.4. Effects of Noise and Bias in Model Parameters 

In practice, the estimates of parameters such as the 
source and background intensity are random variates, 
and are therefore subject to both noise and bias; we now 
derive the consequences of such errors. These effects, in 
contrast with effects in the previous sections, do not al- 
ter the theoretical PDF of intensity. However, they lead 
to similar changes in the model PDF, which mimic the 
effects of an extended emission region. They can thus 
lead to errors in the inferred parameters (7 S and Ard) 
that are distinct from the fundamental limits posed by 
sampling. 

We have seen that the PDF modifications that arise 
from both an extended emission region and decorrela- 
tion of the scintillation within the averaging time are 
naturally represented in terms of partial derivatives of 
P(I;N\G) with respect to G, prior to integration with 
P[G). To express all the results of this section in the 
same way, we again analyze departures using the i.i.d. 
approximation for P(I; N\G). We also give the expected 
bias to 7 S that arises from each type of parameter error. 

4.4.1. Effects of Biased Background or Source Amplitude 

Many effects can bias the estimates of the background 
or source amplitudes. For example, leaked pulse power 
from quadrature downconversion is frequenc y- reversed 
and can be dispersed into the off-pulse region (|Demorestl 
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2007). Moreover, the analog-to-digital conversion of the 
signal requires quantization of the signal; this process 
introduces an offset of the background no i se that varies 
with signal intensity (jJenet fc Anderson! 119981 iGwinnl 
120011 . 

A bias 5I n in the estimated background noise will re- 
sult in a bias of the estimated pulse amplitudes 5Aj = 
—5I n /I s , leading to a bias of the estimated mean inten- 
sity within a scintillation element 5(1 = AI S G + I n ) = 
(1 — G) 5I n . To obtain the incurred PDF modification, 
we apply the i.i.d. approximation and use that P(I; N\G) 
then only depends on the scale I and the degree of av- 
eraging TV. We first expand to leading order in 51, then 
re-express the derivative using AI s d/dI = d/dG, and 
finally integrate G with its exponential weight to obtain 



P(I; N,SI n =0)-P(I;N,5I n ) 



6In 

Ah 



dG (l-G)e 



(21) 

_ G dP(I; N,SI D =0\G) 
dG 



Comparison with Eq. [T7] shows that this incurred mod- 
ification is identical to a change in emission size 7 S — > 
7 S — <5/ n /(2^4is)- A bias in pulse amplitudes is equiva- 
lent, with the substitution <5/ n = — I s 5Aj. Thus, at lead- 
ing order, emission size effects are indistinguishable from 
a biased estimate of the source or background intensities, 
regardless of averaging N. However, because the source 
emission size and the noise bias divided by signal-to-noise 
likely vary differently with pulse amplitude, segregating 
the pulses by strength provides a means to resolve this 
degeneracy. 

4.4.2. Effects of Background or Source Parameter Noise 

We similarly analyze the consequence of unbiased noise 
in the estimated source and background model param- 
eters. This type of parameter error introduces spuri- 
ous modulation in the corresponding model PDF, which 
must be compensated by artificial size effects. We there- 
fore expect a positive bias on the inferred size. 

Such parameter noise is frequently dominated by the 
source power estimate because of the limited sampling of 
the scintillation ensemble in each spectrum. This noise 
5 A will be effectively independent of N, because we have 
assumed that the scintillation pattern remains approxi- 
mately constant during the averaging, so each of the N 
consecutive pulses will be similarly biased. The mean 
intensity within each scintillation element is therefore 
biased: 51 = (G — 1) I S 5A. However, because this pa- 
rameter noise is generally unbiased, the leading order 
correction after an ensemble average is quadratic in 51. 
Again translating to derivatives with respect to G and 
averaging over an ensemble of pulse amplitudes, we ob- 
tain 



P(I; N, 5A=0) - P(7; N, 5 A) 

'5_A^ 2 
A 



(22) 



2 G d*P(I;N,5A=Q\G) 
dGGe dG* ' 



This departure is similar to, though not strictly degener- 
ate with, the effects of a finite emission size. This noise 
5 A increases the inferred dimensionless emission size by 

7s —> 7s + (^(SA/A) 2 ^ /4< The normalized variance is 



given by the reciprocal of the number of averaged scin- 
tillation elements, iVscint ~ B/Az^d, so 7 S experiences a 
positive bias ~l/(4A r sc i nt ). 

The analogous emission size bias from noise in the es- 
timated background intensity is ~l/(4S' 2 Afb), where S 
is the gated signal-to-noise, and Nb is the number of 
samples averaged to estimate the background intensity. 
Observe that parameter noise for the background inten- 
sity, unlike parameter noise arising from scintillation, is 
independent for each spectrum and can therefore be mit- 
igated by averaging over series of pulses. 

Similarly with effects of bias in the model parameters, 
the near-degeneracy with effects of finite emission size 
is broken by comparing multiple analyses. For exam- 
ple, varying the bandwidth used to estimate the model 
parameters changes their noise; the evolution of mea- 
sured residual structure with analyzed bandwidth thus 
precisely quantifies the influence of parameter noise. 

4.5. Effects of Instrumental Distortion 

Many instrumental distortions depend on the signal 
amplitude; they can therefore be modeled by imposing 
appropriate distortions on the distribution of scintillation 
gain P(G). The influence on the PDF can be modeled by 
combining the modified P{G) with the conditional dis- 
tributions P(I;N\G) derived above. Alternatively, the 
modulation index m (Eq. I15p provides a simple measure 
to quantify changes in P(G); we therefore use m to char- 
acterize effects in this section. 

Instrumental distortions can arise, for example, from 
saturation of the observing system. Because strong pul- 
sars can be highly variable and dominate backgrounds 
when they are on, optimal quantization levels are often 
difficult to determine, and even systems with many bits 
may saturate. A simplified model of saturation imposes 
a voltage cutoff on the observed scalar electric field in 
the time-domain: 



x sign(xi) 



if \xi\ < cr D 
else. 



(23) 



The effect of this cutoff is two-fold: power from the 
strongest samples is spread across all channels in the 
spectral domain, and the observed signal intensity is un- 
derestimated. At leading order, the second effect pre- 
dominates, so that the estimated off-pulse noise overes- 
timates the effective on-pulse noise. In terms of the dis- 
cussion in i)4.4.11 saturation leads to a bias 5I n < and 
thus, an inferred emission size that is upward biased. 

Likewise, if the pulse is heavily dispersed, then the sat- 
uration predominantly affects samples with the highest 
scintillation gain G. As a simplified model of this effect, 
we replace the distribution of scintillation gain P(G) = 



e by a truncated version: P(G) 



= e G ( 



(G n 



G) 



K 5 (G — G max ), where 6(x) is the Heaviside step 



function. The corresponding modulation index is 



(G)p 



1 



sinh(G max ) - G 
cosh(G max ) - 1 



< 1. 



(24) 



Comparison with Eq. [T5] shows that even a relatively 
large cutoff G max = 5 gives suppression equivalent to 
modest emission size effects, 7 S ~ 0.014. Furthermore, 
the exponential character of the incurred suppression 
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leads to an extremely rapid onset of saturation artifacts 
once the largest intensities clear the threshold for satura- 
tion; the effects on the strongest pulses may therefore be 
quite pronounced relative to those on the weaker pulses. 

4.6. Expected Resolving Power 

Poisson noise presents a fundamental limit on the re- 
solving power afforded by scintillation, as a fraction of 
the magnified diffractive scale. This limit depends on 
both the gated signal-to- noise ratio S = (A)I s /I n and 
the observational parameters B and T b s . Figure [4] il- 
lustrates the variation of PDF modification with S, for 
N = 1. The modification from extended emission size is 
maximal at zero intensity, with fractional strength 



P(I=0; N=l, 7s ) - P(I=0; N=l, 7s =0) 



P(7=0;iV=l,7 s =0) 

1 e- l ' s 
-5 + (1 + ^EiFI75) 



(25) 



-27s 



where Ei(cc) is the imaginary error function. Equating 
Eq. [25] with the level of Poisson noise r\ in the observed 
PDF then provides an excellent assessment of resolving 
power. 

We apply an approximation for Eq. [25] to obtain the 
following estimate for the achievable resolution, ^/7 s ,min, 
in units of the magnified diffractive scale: 



exp 



0.2 



0.09 (In Sf +0.005 (\nS) 3 



(26) 



This approximation is good to within 10% for 10 
S < 2. Simulations indicate that substituting n 



-3 



< 



— 1.3 



effectively estimates the standard errors from the full dis- 
tribution fits, where iVtot is the total number of sampled 
points used to estimate the PDF of intensity. 

We now apply this result to estimate the resolving 
power for typical observational parameters. A (1 hour)x 
(4 MHz) observation of a pulsar with 5% duty cycle gives 
~10 9 on-pulse samples of the flux density. A gated SNR 
5=1 then gives a 3u detection limit of about 1% of the 
magnified diffractive scale, S = 0.1 gives about 10% of 
the magnified diffractive scale, and S — 0.01 gives about 
50% of the magnified diffractive scale. Since the detec- 
tion limit is much more sensitive to S than to the number 
of sampled points, an analysis of a subset of strong pulses 
may provide optimal resolution. 

5. MOMENTS OF THE DISTRIBUTIONS: THE 
MODULATION INDEX AND SELF-NOISE 

Moments of the intensity distribution provide simple 
diagnostics for scattering and source inference. Fre- 
quently, moments provide optimal estimators of distri- 
bution parameters and are therefore both powerful and 
elegant tools for characterizing distributions. However, 
a blind calculation of moments may be blind to effects 
of RFI, artifacts of finite scintillation averages, or defi- 
ciencies of models, whereas a comparison with the full 
distribution functions might identify such effects. 

We focus on two applications of moments of the distri- 
bution: the intensity modulation index and a practical 
technique for estimating self-noise. For the modulation 
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Figure 4. Intensity PDF and the corresponding modification for 
emission size effects, for N = 1 and varying gated signal-to-noise 
(5). In each case, (A)I B + I n = 2.0. Observe the the plotted 
modification is scaled by both the dimensionless size parameter 
(7s) and S. 

index, we include all sources of noise, as well as evolution 
of the scintillation pattern and spatially-extended source 
emission. For the self-noise estimator, we constrain am- 
plitude variations of the pulsar and assume that there 
is no decorrelation of the scintillation pattern within the 
averaging time (N <C Ard); we then assess departures 
from these assumptions qualitatively. 

5.1. The Modulation Index 

In > j4.2l we demonstrated that the size of the emis- 
sion region affects the modulation index m. However, 
self-noise, temporal decorrelation, and background noise 
also affect the modulation index. These effects are read- 
ily evaluated, cither by analysis of the random variables 
that make up each sample of intensity (Eq. [1]), or by ap- 
propriate derivatives of the corresponding characteristic 
function. 

Here, we analyze only the case for which the pulse am- 
plitudes are i.i.d.: (AjA,-) = (A) 2 + (SA 2 ). This con- 
dition is different than the assumptions for the i.i.d. ap- 
proximation derived above, which assumes the stronger 
condition that the spectral statistics are i.i.d.. Although 
many pulsars appear to have i.i.d. pulse amplitudes, oth- 
ers display striking departures. For example, phenom- 
ena such as nulling and mode-changing break statistical 
isotropy, and the bright est pulses from the Vel a pulsar 
tend to occur in groups (jPalfrevman et al.ll2011[ ). 

To explicitly account for the temporal decorrelation of 
the scintillation pattern, we again assign the autocorre- 
lation appropriate for square-law scattering, quantified 
by the decorrelation timescale (in pulse periods) Ard- 
We incorporate effects of an extended emission region in 
terms of the dimensionless size parameter 7 S , defined by 
Eq.rrj] with the assumption that 7 S < 1. Modification of 
the scattering assumptions will introduce minor scaling 
alterations for these contributions. 
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To leading order in the effects of emission size and tem- 
poral decorrelation, the modulation index is then given 

by 



m 2 = 



_ ((I-In) 2 ) 



- 1 



(27) 



Emission Size Self-Noise Signal-to-Noise 



l-4 7s 



1 \ 2 
1+ n) + N 



1 

~NS 



2 + l 



, (5 A 2 ) , s N 2 

+ 4(1-2^)^^ - (l-4 7s ) 



1 



N(Ay 



6Ar 2 



Intrinsic Variability 



Decorrelation 



Here, we have grouped the contributions to the modula- 
tion index, to explicitly illustrate the relative influence of 
the various effects. We also have used 7s, l = 7 S ,2 = 7s be- 
cause of the degeneracy at linear order. We see that the 
relative contributions of an extended emission region and 
decorrelation of the scintillation pattern give the previ- 
ously derived relationship between their respective PDF 
modifications (Eq. I4.3[) . 

The modulation index corresponds to the zero time 
and frequency lags in the autocorrelation function of in- 
tensity; analysis of non-zero lags eliminates many of the 
terms in Eq. [37] For example, non-zero spectral lags 
eliminate the 'Signal-to-Noise' contribution and mitigate 
the 'Self-Noise' contribution (although intrinsic modula- 
tion induces correlations in the spectral self-noise that 
preserve the contribution). Non-zero temporal lags fully 
eliminate all terms other than 'Emission Size' and 'Decor- 
relation,' which is enhanced. 

5.2. An Estimate of Self-Noise 

The self- noise of a signal is a useful diagnostic, 
which can be used to identify intrinsic variability 
on timescales shorter than the spectral accumulation 
time. Intrinsic modulation induces correlations among 
frequencies, without m odifying the mean spectrum 
(|Gwinn fe; Johnsonll2011[) . If spectra are formed accord- 
ing to the criteria in £12.31 then the power in each spectral 
channel is an exponential random variable. Because such 
variables are completely characterized by their mean, the 
single-channel statistics are therefore immune to intrinsic 
modulation. However, the intrinsic modulation does in- 
troduce correlations in the noise of nearby channels. We 
now derive the expected self-noise in the spectral domain, 
thereby enabling detection of such correlations. 

One method for estimating self-noise is to compare 
pairs of nearby samples that are assumed to be within 
the same scintillation element. We will assume that 
the pair consists of samples with uncorrelated self-noise 
(e.g. pairs from different pulses, or from the same pulse 
with negligible intrinsic modulation). Because self- noise 
is heteroscedastic (see §3.2.3^ . we calculate the noise 

SI 2 = (Ii — I2) 2 /2 in pairs of samples 1%) as a func- 
tion of their mean I = (I1+I2) /2\ we denote the resulting 
noise estimate SI 2 (I). The ensemble average is given by 



J 21 2(1 - I') 2 P(I\ 21 - V- N)dl' 



SI 2 (I) 



J 21 p(r, 21 - 1 1 -, N)dr 



(28) 



where P(Ii, I2] N) represents the probability of sampling 
the pair {I\ , I2}, including an ensemble average over scin- 
tillation. Because we assume that the samples are inde- 
pendent but are drawn from the same scintillation ele- 
ment, P(h,I 2 ;N) = /°° dG P(I%; N\G)P(I 2 ; N\G). 

As usual, we apply the i.i.d. approximation. Although 
the scintillation integral cannot be evaluated in closed 
form, Eq. IA2I readily yields the identity 

,21 

/ 2(I-I') 2 P(I',2I-I';N)dI' (29) 
Jo 

I 2 f 21 
= W~T P{I',2I-l'-N)dI'. 
^ + 2 Jo 

This relationship immediately unveils the simple rela- 
tionship between the estimated mean signal and noise: 



6I 2 (I) = 



N- 



(30) 



In fact, because the identity given by Eq. [29] holds re- 
gardless of the integration over G, this property of the 
estimated self-noise is independent of the degree of sam- 
pling and the properties of scintillation. We see that 
this method for estimating the self-noise gives the same 
form as the exact expression (see £|3.2.3[) . except that 
iV^iV + 1/2. 

If the averaged intensities are from different pulses, 
then pulse-to-pulse variations contribute additional 
noise. If the averaged intensities are from the same 
pulses, then intrinsic variations on timescales shorter 
than i acc induce correlations in self-noise, and thereby de- 
crease the measured noise; if TV > 1 then pulse-to-pulse 
variations within the averaging will increase the noise. 
Each of these effects can be estimated by applying the 
appropriat e substitution TV — » N e ff. Similar tests were 
applied by iGwinn et al.l (|201 If ) to infer short-timescale 
variability (< 300/is) of PSR B0834+06. 

6. SUMMARY 

We have presented observable quantities that describe 
the flux density of a source exhibiting strong diffractive 
scintillation. These observables include both distribution 
functions and bulk indicators, such as the modulation in- 
dex and estimates of self-noise. Our results encompass 
a broad range of physical and instrumental effects, such 
as non-stationary background noise, arbitrary temporal 
averaging, pulse-to-pulse variability, decorrelation of the 
scintillation pattern within the averaging, and the possi- 
bility of spatially-extended source emission, making them 
well-suited for direct comparison with observations. The 
primary benefits of such comparisons include the follow- 
ing: 

• We can completely decouple effects arising from 
decorrelation of the scintillation pattern from those 
of extended emission size by analyzing Nyquist- 
limited data without averaging. 

• Our results enable estimation of the emission region 
sizes of individual pulses and of arbitrary subclasses 
of pulses. 

• For point-source emission, our statistical descrip- 
tion does not require specification of either the na- 
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ture or geometry of the scattering material. Resid- 
uals relative to the point-source model of ^there- 
fore provide robust indicators of either intrinsic 
emission effects (e.g. finite size of the emission re- 
gion), extrinsic effects (e.g. evolution of the scin- 
tillation pattern), or instrumental limitations (e.g. 
finite observing bandwidth). These competing ef- 
fects can be distinguished by varying the degree of 
averaging, the analyzed bandwidth or the chosen 
subset of pulses. 



Our technique does not appeal to extraordi- 
nary scattering geometry (such as is inferred for 
parabolic arcs) or extreme scattering events to in- 
fer the size of the emission region. Furthermore, 
our analysis is effective on short observations (< 1 
hour) and is easily reproducible. The achievable 
resolution limits (see £)4.6[) can be a small fraction 
of the magnified diffractive scale. 



In particular, our results are ideal for describing the 
statistics of low-frequency observations of scintillating 
pulsars. In such cases, the scintillation bandwidth can 
be narrow relative to the observing bandwidth, allow- 
ing an estimate of the scintillation-averaged source in- 
tensity for each pulse. The tight coupling of scintillation 
and source variations is then described by our models, 
which fully incorporate the intrinsic source amplitudes 
and the background noise estimates via on-pulse and off- 
pulse averages. As we report elsewhere, we have found 
excellent agreement between these theoretical expecta- 
tions and observation s of the Vela pulsar at 800 MHz 
(|Johnson et al.l l2012af) . We thereby achieved a spatial 
resolution of ~1% of the magnified diffractive scale - 
about 4 km at the pulsar. Application to other pul- 
sars will enrich the empirical description of pulsar radio 
emission regions, potentially including variation among 
individual pulses and pulse subclasses. 

We thank the U.S. National Science Foundation for 
financial support for this work (AST-1008865). 



APPENDIX 

A. MATHEMATICAL RESULTS 
A.l. Average of Exponential Random Variables with Different Scales 
Consider an exponential random variable /: P(I;I) = I~ 1 e~ 1 / 1 for 7 > 0, where I = (I). The corresponding 
characteristic function is <p(k) = (l — ilk) . The distribution P (/; {Ij}) of the average of N such variables from 
independent distributions with different scales {Ij} is then easily derived by inverting the partial fraction decomposition 
of the product of characteristic functions: 



P(l;{I j })=N 



N 

E 



TN—2 



N 

n 



I>0. 



(Al) 



In other words, the distribution of the average of exponential random variables is a weighted sum of the marginal 
distributions. This sum is the hypoexponential (or generalized Erlang) distribution. 

The representation of Eq. IA1I is singular if any of the scales are equal. In particular, if the random variables are 
i.i.d., each with scale /, then their average is drawn from an Erlang distribution, which is simply a gamma distribution 
with integer shape parameter: 



P(l;lN) 



N 



N 



I 



JV-1 



J>0. 



(iV-1)! I N 

iGoodmanl (|1985l ) discusses both these distributions in the similar context of polarized thermal light. 



(A2) 



A. 2. Distribution of Correlated Exponential Random Variables 

We now approximate the distribution of N correlated exponential random variables {h}. The exact distribution 
follows by first constructing a complex multivariate normal distribution with the desired covariance matrix; the ex- 
ponential random variables are the squared norms of each circular complex Gaussian random variable. We restrict 
ourselves t o the case in which each random variable has unit mean. The characteristic function of P ({h}) can then 
be written (|Mallildl2003T ) 



1 



det [I-idiag{fcJS]' 



(A3) 



Here, [E] 4j = ^(IJ,) - 1 = 
We now expand in e = J 



where J is the NxN unit matrix: [J]y = 1 V Since [e]ij ~ (1 — Fy) /2, we can 



Intensity Scintillation Statistics 



13 



approximate the characteristic function to leading order as 



v({fci};=) = 



(A4) 



The inverse Fourier transform of c/j({fci};S) then gives our approximation to P({7j}). Although this inverse can be 
written in terms of delta functions and derivatives, the most natural form expresses the action as a weight function: 



J d N h f({Ii})P({Ii}) = £ dl /({/})e- 7 "X>- r «) /e ~' 7 + ° (j 1 - r «l 



(A5) 



Here /({/}) denotes /({/,; = J V i}). 



A. 2.1. Example: Bivariate PDF of Correlated Exponential Random Variables 

The joint distribution of a pair of correlated exponential random variables is particularly simple. In this case, the 
characteristic function (|A3I) is easily inverted: 



P(hJ 2 ;T) = - J— 1 l^^ JIJ-A exp 



1 



(i-T) (/l 



/i,/ 2 >o. 



l -r~ u Ui-r) 

Here, T — (I1I2) — 1, and 7n(a;) is the modihed Bessel function. Using Eq. IA51 we can approximate Eg. IA6l as 

2-7 



P(I 1 ,I 2 ;T) = e- 1 U(AI) + (1-T) 



where 7 = (h + 7 2 )/2 and A7 = 7i - 7 2 . 



-<5(A7) +7(T(A7) 



(i-iy 



, 7 > 0, 



(A6) 



(A7) 



A. 2. 2. Convolution of Correlated Exponential Random Variables 

The gamma distribution is an excellent approximation for the sum of independent exponential random variables that 
are not necessarily isotropic (see £|3.2[) and can also be used to approximate the convolution of correlated exponential 
random variables. This approximation is most effective for correlations that are particularly strong or weak. For 
example, consider the average of two exponential random variables, each with unit mean, and with covariance (Iih) = 
1 + r. The gamma distribution that matches the mean and variance of the exact average has a corresponding effective 
number of degrees of freedom N c g = 2/(1 + T). These distributions agree exactly for the limits r — > 0, 1, but 
suffer sizeable discrepancies for intermediate correlations. For example, the third central moment of the average of 
two correlated exponential random variables is (1 + 3r)/2 but is 2/N^ s = (1 + T) 2 /2 for the corresponding gamma 
distribution - a difference of ~10% for intermediate correlation. 

More generally, for N averaged exponential random variables with correlation matrix Tij, as defined above, and unit 

mean, the appropriate number of degrees of freedom based on variance matching is N e g = N ^1 + 2A _1 J2i<j ^ij^j ■ 

Even in the case of large N, the disparity in skewness between the exact average and the appropriate gamma distribution 
remains. For instance, if = T for i ^ j then, 



fJ-3 



2r 3 / 2 + f (l-Vf 



2 

N 2 



2T 2 



4T 

N 



(i-r) 



2 

TP 



(i-r) 



1 -3r 

2 



2r 3/2 ) Exact Average, 

Approximation with A c ff = 1+ ^ N 



N 



(A8) 



-i)r' 



where ^ denotes the third central moment. Again, we see that these expressions agree for T 
differ substantially for intermediate correlation, even as N 00. 



0, 1 as expected, but 
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